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Abstract 

We present a model for semiflexible polymers in Hamiltonian for- 
mulation which interpolates between a Rouse chain and worm-like 
chain. Both models are realized as limits for the parameters. The 
model parameters can also be chosen to match the experimental force- 
extension curve for double-stranded DNA. Near the ground state of 
the Hamiltonian, the eigenvalues for the longitudinal (stretching) and 
the transversal (bending) modes of a chain with N springs, indexed by 
p, scale as Aj, ~ (p/N)'^ and A* ~ p^(p — 1)^/A^ respectively for small 
p. We also show that the associated decay times Tp ~ {N/pY will not 
be observed if they exceed the orientational time scale ~ for 
an equally-long rigid rod, as the driven decay is then washed out by 
diffusive motion. 

1 Introduction 

There are many models for the equihbrium behavior of long polymers but few 
for their dynamics. The celebrated Rouse model |T] describes the temporal 
behavior of the chain by a set of (Rouse) modes in a harmonic potential. Such 
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a Rouse (or fantom) chain is however too flexible to describe the behavior 
of biological polymers. For instance the persistence length of a Rouse chain 
is zero, while e.g. double-stranded (ds) DNA persists in its orientation over 
some hundred base pairs. It is easy to dress the Rouse model with a term 
that suppresses bending of the chain, while preserving the dynamics in the 
form of a set of independent Rouse modes. However, the behavior of such 
a chain under an external force, i.e. the force extension curve, does still 
not correspond to that of the semifiexible biopolymers, as it responds to the 
force as a harmonic spring, while the semifiexible chain has a characteristic 
rapid increase of the extension for small forces, crossing over to a slow linear 
increase with the force. Clearly Rouse-type chains have too little rigidity 
against longitudinal forces. 

An opposite approach is to represent the biopolymer as a worm-like chain 
[21 13] , which is a continuum description with the appropriate stiffness against 
bending. The worm-like chain (WLC) keeps a fixed contour length and 
only undergoes transversal undulations. A microscopic realization of the 
worm-like chain is the model in which the monomers are connected by rigid 
bonds. The bending stiffness is incorporated in the hamiltonian by a term 
suppressing the angle between two successive bonds. The rigidity of the 
bonds guarantees the invariance of the contour-length under the possible 
motion of the monomers. We refer to this model as the "jointed chain" (in 
analogy with the well known freely-jointed chain p]). A refinement is to let 
the bond length vary with the external force. Such an "extensible jointed 
chain" can very well reproduce the measured force-extension curve of bio- 
polymers by an appropriate choice of the bending stiffness. A model close to 
the extensible jointed chain has recently been introduced by Kierfeld et al. 
[S], who extensively discuss its equilibrium properties. 

However, the dynamics of jointed chains is quite involved, as the motion 
of the monomers is strongly constrained by the rigidity of the bonds, which 
introduces unphysical long-range correlation. E.g. the motion between the 
end monomers of the chain is correlated due to the constant length of the 
contour. 

Thus a polymer model that can treat the dynamics of biological polymer 
fragments, such as occuring in networks, is quite welcome. These fragments 
are of the order of the persistence length, rendering a continuum description 
not adequate and a description on base-pair level necessary. In this paper 
we propose such a model combining to a large extend the mathematical 
simplicity of the Rouse model, while describing accurately the experimental 
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force-extension curve. 

The hamiltonian for our flexible and extensible chain reads: 

\ N N-l 

^ = 2 Ed^nl -df-f^T. ■ u„+i - F ■ L. (1) 

n=l n=l 

Here u„ is the bond vector between monomer n — 1 and n and L the end-to- 
end vector 

N 

Un = Tn — fn-l, L = rjy — Fq = ^ U„. (2) 

n=l 

r„ is the position of the n-th monomer {n = 0,1, ■ ■ ■ , N). F is a force 
tending to orient and stretch the chain. The first term in ([1]) is a single-bond 
harmonic interaction, providing the length scale d and the energy parameter 
A. The second term is a nearest-neighbor-bond interaction giving the model 
a bending stiffness measured by n. The above defined hamiltonian is the 
simplest of a class of hamiltonians with bond interactions, not only between 
nearest neighbors as in ([1]), but also longer-ranged interactions. If we set 
d = the model reduces to that of Marques and Frederickson [BJ (with slightly 
different boundary conditions). As we mentioned above, such a Rouse-like 
model does not have an acceptable force-extension curve. We will argue that 
a finite d changes the properties from being rather unrealistic to realistic. 
The difference between our model and that of Kierfeld et al.[5] is minor 
as far as the equilibrium properties are concerned. However the difference 
is major with respect to the dynamics. While in |S] the bending energy 
only depends on the angle between the bonds, in our case it depends on 
the bond vectors. The angles are dynamically very inconvenient parameters: 
changing one angle while keeping the others fixed, involves the motion of a 
whole segment of the chain. In our model monomer positions can be changed 
independently, only changing the incident bond vectors of that monomer. 

In this paper we discuss the mechanical and statistical properties of the 
model. We formulate the dynamical equations for the monomer positions, as 
well as those for an equivalent description in terms of the modes of the system. 
We also give a few implications of the mode equations. In a subsequent study 
we design an efficient scheme for the solution of the dynamical equations and 
report on extensive simulations of the model. 
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2 The parameters of the model 



By varying the parameters d, X and k we cover a wide range of physical 
systems. We can ehminate the value of a finite d from the problem by scaling 
u„ with d. Then the hamiltonian reads 

Xd^ iV N-l N 

n = — Xl(|un| - 1)^ - K.d^ u„ ■ u„+i -dY u„. (3) 

n=l n=l n=l 

Thus d can be absorbed in the parameters A, k and F. Therefore d does not 
show up directly in the expressions for the static and dynamic properties. To 
translate the properties to a real chain one has to multiply the bond vectors 
with a well chosen d. The form (|3]) shows that c/ = is a singular point of 
the hamiltonian ([T]). Amongst others, the case = has a lowest energy 
state in which all bond lengths vanish. 

As our main interest is in describing biological polymers, which have a 
fairly rigid bond length and a large persistence length, we consider large 
values of A and k. So it is useful to consider the ratio between A and k as 
characteristic for the model. As we will see their common magnitude couples 
to the temperature. So we write 

K = vX (4) 
and the hamiltonian gets the form that we will use in this paper: 

with the abbreviation 

f = F/(Ac/). (6) 
Equilibrium is governed by the ratio T-L/ksT, which we write as 

jL = n_ (7) 

keT T* ^ ' 

where the reduced temperature T* equals 

T* = keT/iXd^). (8) 

Thus the combination Xd'^ is absorbed in the reduced temperature T*. We 
will discuss the model in terms of the dimensionless parameters T*, u and /. 
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Large values of A manifest themselves in low values of T*, implying a fairly 
constant bond length. Assuming a constant bond length reduces the model 
to the jointed chain and indeed the force-extension curve of our model closely 
follows that of the extensible jointed chain. 

For the analysis of the model we split the hamiltonian in three parts 
according to the power in which the bond vectors (or monomer positions) 
occur: 

H* = HI + HI + N/2. (9) 

The first term, containing the quadratic terms in the bond vectors, has the 
form 

N N-l 

'^h = 7^Y.^l-^Y.^n- U„+i. (10) 

It is instructive to compare this hamiltonian with the one used by Marques 
and Frederickson [6], which, using our parameters, reads 

1 - 2z/ ^ u^-^ 

n = — — E «n + IT E - ^n+i)'- (11) 

^ n ^ n=l 

One observes that stability requires the bound u < 1/2, which we also have 
to impose on Close inspection shows that ffTUl) and flTT]) are the same, 
apart from two boundary terms: ul and u% occur in flTUl) with the coeffi- 
cient 1/2, while in (fTT!) they have the weight (1 — z^)/2. For the polymer 
properties this difference is of little importance, but the form ffTOl) permits 
exact diagonalization in closed form, while for (fTTl) one has to rely on numer- 
ical diagonalization. We will use "H^ as expressed in terms of the monomer 
positions 

1 ^ 

Hfi = — Hm,n I'm, ■ r^. (12) 

m,n 

We will show that the matrix ifm,n can be exactly diagonalized. It would be 
the full (force-free) hamiltonian if we were to set d = 0. 

The term Hs represents the stabilizing force and has the form 

N 

K = -EK + f -Un]. (13) 

n=l 

The first term in f|T3l) . involving the contour length of the chain, plays an 
essential role in stabilizing the chain. It prevents however an exact diagonal- 
ization of the model, which is possible for the harmonic terms. The second 
term in f[T^ accounts for the infiuence of the force on the chain. 
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The last term in ([H]) is a trivial constant having no influence, neither on 
the equilibrium nor on the dynamics; we will pay no further attention to it. 



3 Ground state Properties 

The dynamics of the model is our main interest. There are a number of 
equilibrium properties which are important for the dynamical analysis. We 
summarize them here. 



3.1 The ground state configuration 



The ground state or lowest energy conflguration prevails at low T*. With 
a non-vanishing external force, the chain starts to align itself with the force 
and the configuration becomes a straight line. All bond vectors point in the 
direction of f and may be written as u„ = M„f . The values of the length Un 
follow from minimization of the energy, leading to the equations 



/ Ui — VU2 



— UUn-1 + UN 



1 + /, 
1 + /, 

1+f.J 



(14) 



Note that the solution for / 7^ follows from that for / = by multiplying 
all M„ by 1 + /. A set of equations of the type ffT^ can be solved by the 
standard technique of the generating function. The solution reads 



1+/ 

1 - 2u 



cosh{a(A^ + 1 - 2n)} 
cosh{a(A^ + 1)} 



(15) 



with cosh(2Q;) = 2/v. For the bulk value we take the bond n = N/2 in the 
middle. The hyperbolic cosine in the numerator then obtains the argument a, 
while the denominator is exponentially large in A^. Leaving out exponentially 
small contributions the bulk values equals 



Kf) 



1 + / 

1 - 2u' 



(16) 
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The bonds become somewhat shorter near the ends of the chain. The value 
b = 6(0) is an important (numerical) value by which d has to be multiplied 
in order to match the experimental bond length a. In other words d = a/b. 
We find that for small T*, 6 is still close to the value f|T6|) . 



3.2 Eigenmodes near the ground state 

The ground state energy responds harmonically to small deviations, yielding 
a quadratic increase in the energy. We first discuss the force free case. One 
can find the Hessian matrix by twice differentiating the Hamiltonian with 
respect to the positions of the monomers. As the dependence on the monomer 
positions is exclusively through the bond vectors, we use the rule 



d 



d 



d 



5r„. dur. 



du 



n+l 



and define the bond matrix 



(17) 



(18) 



dUmdUn 

The indices in the matrix B run between 1 < m,n < N. We extend the 
matrix by putting the elements equal to zero, whenever one of the indices 
equals or + 1. The Hessian matrix M is then obtained as 



B 



r?i,(n+l) 



B 



(m+l),n 



+ B 



(r?i+l),(n+l)- 



(19) 



Mm,n runs in the interval < m,n < N. 

Let us first consider the part of B which is due to the contribution 'Hh- 
This part is diagonal in the vectorial indices and the scalar matrix multiplying 
the unit tensor I reads 



/ 1 



B„ 



-V 



-V 

1 -V 

-V 1 








-V 



(20) 



Using the property 
min + l)7r 



sm 



Ar + 1 



+ sin 



m{n — l)7r 
Ar + 1 



2 cos 



ivTT 



sm 



mrni 



(21) 
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it is easy to see that the eigenvalue equation has the form 

pniT 



N 
n=l 



sm 



N + 1 



with 6p given by 



On 



2z/cos 



. / pnn 
Ur, sm 



pn 



(22) 



(23) 



Note that fl22ti23p are also satisfied for the ends of the chain {m,n = 0,N) 
(which do not correspond to bond indices), since we can trivially extend the 
summation to the interval < m, n < A^. 

Using Eq. fl22|) and similar trigonometric relations, we arrive at the eigen- 
values for Mm„ 



N 



EM 



pn 



COS 



n=0 



(n + l/2)7r 



(p cos 



(p+l/2)7r 
N+1 



N+1 

with (p given by 

(p = 2 1 — cos 

Next we consider the infiuence of the stabilizing force Since 

Q2 N 



\N + 1 



9 p. 



evaluation of the derivative yields 



B 



I 



u, u. 



(24) 



(25) 



(26) 



(27) 



This tensor projects out the longitudinal eigenmodes and is unity for the 
transversal modes. Longitudinal and transversal relate to the direction of 
the vector Uj. Thus the longitudinal modes are given by f l22ll23|) . We put 
this in formula for the decay constants Xp of the longitudinal mode p 



- Cp- 



{2i 



The two transversal modes are corrected by the term (127|) . which is diag- 
onal in the matrix indices. The diagonal contribution to B of the quadratic 
terms in the Hamiltonian equals 1. So if we define 



(1-lM), 



(29) 
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qi becomes the diagonal matrix element. The transverse modes for B follow 
from the the diagonalization of this matrix, which has to be carried out 
numerically. One gets an impression of the eigenvalue spectrum by using the 
bulk value 

l/ui ~ 1 - 2z/ (30) 

for each diagonal element. Then the modes of the bond matrix have again 
the shape and the eigenvalues are 



1 — cos 



iV+ 1 



(31) 



This formula covers the spectrum closely, provided that we replace for the 
low modes p by p — 1 in the argument of the cosine. Note that upon this 
substitution the eigenvalue for p = 1 vanishes. This is a general property of 
the eigenmode in the force free case, due to invariance of the energy under 
an overall rotation. 

The dynamical matrix then has transversal eigenmodes similar to (!25il 
and the eigenvalue equals 



1 — cos 



N +1 



el. (32) 



An external force changes the shape of the modes and the eigenvalues. 
The transversal mode for p = 1 no longer vanishes. The effects of a force 
have to be incorporated numerically. 

3.3 The persistence length 

The persistence length derives from the orientation correlation function 

Oij = (ui-Uj). (33) 

Asymptotically Ojj decays exponentially and the exponent defines the per- 
sistence length. It is dominantly dependent on v (or k). For the large values 
of A (or small values of T*), one gets a good impression of the orientation 
correlation function, by taking the bond lengths as rigid, i.e. by studying the 
extensible jointed chain, for which the orientation correlation function can 
be calculated exactly (see Appendix One finds for the persistence length 
Ip 

^ = -— — — - with r = b^u/T*. (34) 

a log[coth(r) - 1/r] ' ^ ' 
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Note that only k determines this persistence length, since A, which is hurried 
both in V and T*, drops out. So a given persistence length sets the ratio 
V jT* . This expression gives the persistence length as the number counting 
the monomers over which the orientation persists. 



3.4 The force-extension curve 

The most informative quantity is the force-extension curve. It is the response 
of the end-to-end vector L on the external force. It follows from the partition 
function 

= Trexp-7{/(A;6T) (35) 

as 

(L) = (Eu„> = %^. (36) 

n=l '-'^ 

This quantity is difficult to simulate, since it fluctuates strongly with a long 
relaxation time. The limits of weak and strong forces are however easy to un- 
derstand. For weak forces we may treat the external force as a perturbation, 
with the lowest order result 

logZ^logZo + ^(([f-L]V (37) 

The subscript indicates the forceless values. Thus the response is linear in 
f: 

(L) ~ (LL)o ■ f/T* = f (L2)o/3T*. (38) 

In the forceless average only the diagonal components survive. The mean- 
squared- average of the end-to-end vector can be related to the persistence 
length for chains many times longer than the persistence length: 

(L2)o = 2iV6(0)/p. (39) 

The combination = Nb{0) is the contour length and thus one gets for the 
initial slope in the large N limit 

M ^ 3kt. (40) 

Le ST* ^ ' 

One sees from this expression that for low T* the slope of the rise is steep 
for two reasons: the small factor T* in the denominator and large Ip in the 
numerator. 
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In the large force limit we may use the result (fTB|) . The chain will be 
stretched in the force direction and the bond will be elongated accordingly, 
with the result 

S^^(l + /)f. (41) 

As a note in passing we point out on the basis of equation f|T6l) that the 
ground state gives this result for all forces. Thus for f = there is a residual 
value for the extension. This holds only for T = as for all T 7^ the 
end-to-end vector will vanish for / = 0, due to rotational invariance. Indeed 
the point (/ = 0, T = 0) is singular and the value of (L) depends on the way 
this point is approached. 

The complete curve will be a smooth transition between the fast initial 
rise (till L is of the order Nh) and the slow rise of (1221) . manifesting itself 
when / becomes of the order 1. Remember that / is the real force F divided 
by \d. Thus forces / of order 1 are huge on a real scale when the polymer is 
stiff or A large. In between, there is a quasi-plateau with the value Nh. We 
have shown the force-extension curve for our model in Fig. [T]for a number of 
chain lengths A^. 



3.5 Choice of the parameters 



The model parameters v and T* can be chosen such that the force-extension 
curve of the model is in good agreement with that of biopolymers. Wang et 
al. [7] give the following empirical formula for biopolymers, such as dsDNA: 

-2 



Flp 



1 



F 



1 

-4 + 



{1^ 



F 
K'o' 



(42) 



This equation contains two polymer-specific parameters: the persistence 
length Ip and the constant Kq. As we have mentioned earlier, the exper- 
imental Ip has to be combined with the bond length a to a dimenionless 
quantity for which our model gives the value For stiff bio-polymers T* 
is small, making r large and Ip/a equal to r. So as first equation we get 



6V 



(43) 



The other experimental quantity Kq is turned into a dimensionless combina- 
tion z: 



11 



This constant manifests itself only for large forces. The behavior of the 
extension for large forces is somewhat hidden in the implicit equation (142 p . 
In order to make it explicit, we write the equation in parametric form. We 
introduce the auxiliary parameter s and for the force parameter x: 

Thus X is a function of s alone according to equation ( H2i) . On the other 
hand the definition of x relates it to F/ Kq as 

F . . {L) X 

x=-:—rz, such that ——= s -\ . (46) 

Kq Lc rz 

The parameter s runs in the interval < s < 1; small s corresponding to 
small forces and s close to 1, to large forces. The last equation (H6|) gives the 
extension as function of the parameter s, while the last equation ( HSj) gives 
the force as function of s. 

For s = 1 — e, with e small, one has, according to ( l45l) x ~ l/(4e^) such 
that X ^ s. Then the force-extension curve obtains the form 

^=.1 + -. (47) 
Lc rz 

We compare this with the strong force behavior of our model as given by 
( H2|) . Relating the experimental force to our parameter / we have 

X(P T* aT* T*' ^ ' 

The two expressions (1471) and (HT!) are the same if 

^ = ^- (49) 

We now have two equations, and (H^ . for the two hamiltonian pa- 
rameters T* and u. Solving them leads to the explicit relations 

T* = '^L±1 and u = ^—. (50) 
z^ 2r + z ^ ^ 

For dsDNA the persistence length is about 40 nm. Wang et al. ff] report 
for the force constant a value around Kq = 1200 pN. So with a monomer 
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distance a = 0.33 nm and at room temperature /c^T = 4pNnm, we get as 
representative values 

r = 120 and z = 100. (51) 

This gives for u and T* as model parameters for dsDNA 

T* = 0.034 and u = 0.35. (52) 

In Fig. [T]we show that simulations of our model for large N give a force- 
extension curves in excellent agreement with the curves of the extensible 
jointed chain for finite chain-lengths N = 64, 128 and 256. The extensible 
jointed chain can easily be extrapolated to = oo and this curve follows 
closely the empirical curve describing the force-extension of bio-polymers. So 
our model is capable to describe bio-polymers that follow the empirical curve 
of HI- 



4 The dynamical Equations 

We treat the dynamics of the polymer in the so-called large viscosity limit, 
where effects of dynamical inertia of the monomers are absent and the sol- 
vent is represented by a noise source uncorrelated in space and time. So 
correlations in the solvent motion due to hydrodynamical motion are omit- 
ted. To incorporate these correlations in the noise is another project. Then 
the dynamical equations for the monomers are the Langevin equations for 
particles in a high viscous medium: 

(ir„ 1 &H , . , . 

Here ^ is the friction coefficient and g„ is a random noise force, with corre- 
lation function 

(gn(t) gm(t')) = ^ I<^n,m5(t - t'), (54) 

where I is the unit tensor. The strength of the correlation is determined by 
the fluctuation-dissipation relation. 

We make the equations dimensionless by introducing a reduced time r: 

r = tA/e (55) 
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The equations then obtain the form 
dvm ^ 

—p- = - XI Hm,nrn{r) + h„(r) + g„(r) + f(5„,o - ^N.n), (56) 
hm derives from the stabihzing term Tig, given in f lT^ . reading 

— 1 o ~ l^ml ~ |U,n+l|- J 

u„ is the unit vector in the direction of the bond n. The dimensionless 
random forces have the correlations 

{gn{r)gUr'))=2T*I6n,mS{r-r'). (58) 

The center-of-mass position, defined as 

1 ^ 



N+ . n 



Y E rn, (59) 



is not affected by the interaction between the beads and is only influenced 
by the random forces. It fluctuates due to the random forces as a Brownian 
particle. 



5 Rouse Modes 

Diagonalizing the matrix Hm,n gives modes relevant for the system. Since 
H is symmetric and positive definite, all its eigenvalues are real and non- 
negative. This means that all the modes, but one, are associated with a 
decay exponent (p. Only the mode corresponding to the center-of-mass has 
a zero eigenvalue and decouples from the dynamic equations (the sum over 
the columns of -ffm.n vanishes). 

The eigenmodes Rp of the matrix H induce an orthogonal transformation 
of the positions r„. The label p runs through the integers p = 0,1, . . . N. For 
p = we have a mode corresponding to the center-of-mass 

Ro = (iV + 1)1/2 R,^. (go) 
The other modes are given by the expression 
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The modes are normalized such that the transformation from r„ to Hp is 
orthonormal. For that reason we have put the factor in front of Hem in (160 p . 
One reckognizes the modes Rp as the modes of the Rouse model. The Rouse 
modes have generally turned out to be an important tool for analyzing the 
properties of a polymer chain [8l |9]. The Rouse modes diagonalize the part 
HI of the hamiltonian. The eigenvalues are, acoording to fl2^ . 



COS 



pn 



N+1 



2z/cos 



/ pn 



\N+1 



(62) 



showing that the slow modes vanish as (p/N)"^ for small p. 

As the Rouse modes diagonalize the dynamic matrix H, it can be written 

as 

N N 

Hm,n rm ■ Tn = YCpRl, (63) 
m,n=0 p=0 

with (p given by ( 162|) . We point out that Rp is an exact eigenmode of Hm,n 
because of the special structure of the boundaries of the the matrix. 
The inverse transformation reads 



Rem ~l~ 



1 



1/2 N 

p=l 



{n + l/2)p7r \ 



(64) 



The advantage of these linear relations between the Rouse modes and the 
monomer positions, is that we can switch between either of two representa- 
tions of the configuration to whichever is more convenient. 

By applying the transformation to the equations fl55|) they become 



dRp 

dr 



-CpRp(r) 



Gp{r). 



(65) 



Here Hp is the transform of h„ defined in (I5 



1/2 N 
n=0 



[n 



+ l/2)pvr\ 



The external force on the Rouse modes gets the form 



(66) 



pn 



f. 



(67) 
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As the transformation to Rouse variables is orthogonal, the correlation 
function of Gp equals 

(Gp(r)G,(r')) = 2T*I5p,,5(r-r'). (68) 

Without Hp the Rouse modes are exact eigenmodes with decay coefficient 

We make the scheme more explicit by inserting fl571) into fl66|) and com- 
bining the two terms, yielding 

N 

Hp(r) = -E^P."^"M' (69) 



n=l 



where the matrix Mp „ is defined as 



1/2 



The bond vectors are likewise expressed in terms of the Rouse variables using 
(El: 

N 

u„ = -5]R,M,,„. (71) 

q=l 

By these steps the relation between the stabilizing force and the Rouse modes 
has been made operational. We note that the operations have the character 
of a discrete sine transform and therefore the technique of the fast fourier 
transform can be applied in both relations fl69l) and flTTl) . 

In concluding the discussion of the Rouse modes, we eleborate on the 
relation between the modes of the bond matrix Bm,n, introduced in fl20|) and 
the Rouse modes. To make this relation more transparant, we introduce 
modified Rouse modes as 

Rp = -2 sin (—^^] Rp. (72) 



.2(A^+1), 



Then relation flTTl) can be written as 



iv + i; ^1 \N + i 
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The transformation from the Rp to the u„ is also orthogonal and self-dual: 
R 



up 



—i) 5-(]vtt)"- p*' 



p=l 



This provides actually a short-cut to derive (1621) . The harmonic energy is on 
the one hand with fl20l) given by 

N ^ _ 

Bm,n Urra " U„ = ^^p Rp ■ Rp ("^5) 

and on the other hand it also equals 

N N 

^ ^ Hm,n ' I"n ^ ^ Cp -fi-p ' Rp- l"^^) 

m,n p=l 

With relation fl72|) between the modified Rouse modes and the Rouse modes, 
relation (|62|) immediately follows. 



6 Rotational diffusion 

Apart from the Rouse modes we encountered in section 13.21 the mechanical 
modes of the harmonic approximation near the ground state. The ground 
state is relevant for polymers short with respect to the persistence length. 
Then the majority of the configurations will be similar to a ground state. 
The ground states are degenerate with respect to their orientation. A small 
force lifts this degeneracy. 

Both the Rouse modes and the mechanical modes are linear orthogonal 
transformations of the monomer positions and can therefore be expressed 
in terms of each other. The modes are similar; in fact the longitudinal 
mechanical modes are Rouse modes. Also the Rouse mode p is dominantly 
present in the representation of the transversal mechanical mode p in terms 
of Rouse modes. Nevertheless their role is quite different. Each transversal 
mechanical mode is an exact eigen-mode of the system with an eigenvalue 
quite different from the decay constant of the corresponding Rouse mode. 
Since mechanical modes refer to small deviations from the ground state, it is 
clear that they do not form an adequate representation, if the ground state is 
not dominant (as is the case for polymers longer than the persistence length). 
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The orientation of the ground state is best represented by the end-to-end 
vector L defined in (15^ . It can be expressed in terms of Rouse modes as 



1/2 



E 

p=l,3,- 



COS 




(77) 



showing that the low p modes are most important for the orientation. Like- 
wise the mechanical modes with low p will contribute dominantly in the 
expression of L in terms of the mechanical modes. 

The dynamical equations can also be formulated in terms of the me- 
chanical mode- amplitudes Ap(r). Here the superscript a distinguishes the 
longitudinal (a = /) and the transversal (a = t) modes. As the modes are 
exact eigen-modes, the equations read 



dr 



-x;A;iT) + G;iT) 



(78) 



with A's from section [221 Note that the equations of the mechanical modes 
for different p are independent, since the noise forces G'p(r) are again un- 
correlated. Eq. (178|) is a standard Langevin equation for a variable with a 
systematic force Xp and a random force Gp. The corresponding Fokker-Planck 
equation for the distribution is explicitly soluble [TU]. The decay constant 
determines the time scale on which the memory of the initial value decays. 
For longitudinal eigenvalues we have an exact expression, following from 
(1251) and (El as 



a; = 2 



1 — COS 



piT 



N+1 



1 — 2v cos 



/ p-K 



The slow-mode eigenvalues Aj, decrease with as 



A' ~ 

p 



2u 



P^Tl'^ 



(79) 



(80) 



For the transverse eigenvalues we have no closed expression, but they are 
well approximated by ( l32l) and ( l3Tl) . 



1 — cos 



/ p-K 



\N +1 



which gives the slow modes the form 

p^ip 



cos 



N +1 



(81) 



(82) 
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In Fig. [2]we show their behavior as function of mode index p and length A^. 
The striking feature is the large difference in decay coefficient for the slow 
longitudinal and transversal modes. The decay as {N/p)~^, characteristic 
for the worm-like chain, has often lead to the conclusion that the slowest 
modes in the system show correlations over a time scale increasing as {N/p)^. 
This is however not the case, a fortiori not for transversal mode p = 1, 
having a vanishing decay constant. If one starts from a ground state with a 
certain orientation and looks at the evolution on a time scale of order A^^, 
then the longitudinal modes have come to equilibrium, while the transversal 
components still fluctuate away from the chosen ground-state configuration. 
In fact they cause a diffusion in orientation, while keeping the shape of the 
polymer more or less fixed. The problem of diffusion of a rigid rod has already 
been considered by J. M. Burgers in 1938 [Hj. He calculated the rotational 
diffusion constant Dj. for a rigid rod with moment of inertia / as 

where we have inserted the moment of inertia of the ground state / = 
a^A^^/12; here, we recollect that a is the distance between monomers. This 
corresponds with a rotational relaxation time 

2Dr^ 24T*' ^ ^ 

which is indeed of the order A^^. Thus on the scale of r^, the orientation of 
the original ground state has diffused away and the slow time scale 1/A* of 
the transversal modes will not be observed. 

It is interesting to compare the decay times of the transversal modes 
with this rotational- diffusion time scale Tr. Using f l5^ we find 

1 _ fNa\ 



24 ^ « 1> P> 1- (85) 



This formula seems to indicate that is a factor larger than Tr. However, 
relation (l85ll assumes that the persistence length Ip exceeds the length Na. 
So the autocorrelation time for the amplitude of mode p > 1 is significantly 
5ma//er than x^,. Thus, the scaling of the decay time proportional to {N/pY 
should not be seen as an indication that semifiexible polymers have very long 
times for small p, but as an indication that their high-p degrees of freedom 
decay very quickly. 



19 



7 Discussion 



We have discussed properties of a polymeric chain characterised by a simple 
hamiltonian, in the high viscosity limit with uncorrelated noise. We feel that 
our model fills a gap in the spectrum of models with respect to the dynamics 
of polymers for the following reasons. 

1. With only two adaptable parameters T* and z/, the model interpolates 
between a worm-like chain and the Rouse model. Depending on the 
length of the polymer, short or long with respect to the persistence 
length, the properties are either more worm-like or more Rouse-like. 
This is important for the dynamics of fragments of biological polymers 
such as dsDNA, which neither are worm-like nor are well described by 
a Rouse model. The combination v/T* controls the persistence length 
and with T* — one approaches the worm-like chain. 

2. The empirical force-extension curve for biological polymers as proposed 
by Wang et al. [7] contains also two adaptable parameters: persistence 
length Ip and the force constant Kq. In Section 1X51 we show how to make 
a map between Ip and Kq on the one hand and our model parameters 
T* and v on the other hand. The correspondence between the force- 
extension curve of our model and that of Wang et al. is excellent, as 
Fig. [1] shows. We point out that our model also enables to discuss the 
force-extension of short polymer fragment as occurring in networks. 
The appendix contains the details of such a finite-size analysis. 

3. The equilibrium properties of the model are close to that of the model of 
Kierfeld et al. [S] . The difference is that we use for the bending energy 
the inner product of two adjacent bond vectors, while (commonly) the 
cosine of the angle between the bonds is used. Our model allows to have 
the monomer positions as unconstrained dynamical variables, whereas 
it is awkward to formulate the dynamics in terms of the angles. 

4. The dynamics of the model admits a description in terms of modes, 
showing clearly the wide spectrum of time scales from short for high- 
indexed modes till very long for the important low-indexed modes. Two 
types of modes exist. For polymers much shorter than the persistence 
length, the mechanical modes (described in Section 13. 2p are relevant. 
For longer chains the Rouse modes are the more convenient representa- 
tion. Since the modes have a wide spectrum of properties we have the 
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opportunity to focus on slow and large scale phenomena and to sup- 
press the role of the fast short-range behavior, which normally greatly 
slows down the calculations. 

The model has several options for improvement. One could introduce 
further bond interactions to make the model more suited to described specific 
polymers. In the same spirit one could add hydrodynamic interactions like in 
the Zimm model [121 US] . In a subsequent paper we describe further results 
of simulations of this model for realistic parameters T* and u applying to bio- 
polymers. We there show how the simulation can be speeded up by orders of 
magnitude permitting to study chains of the order of 1000 base pairs, which 
makes the Rouse modes a convenient description of the dynamics. 

Acknowledgement. The authors are indebted to Debabrata Panja for 
numerous stimulating discussions and providing the simulation data of Fig. 



A The extensible jointed chain 

In this appendix we discuss the properties of the chain in which the bonds 
all have the length b{f) = 6(1 + /). The / dependence of this bond length 
derives from the stretching of the bonds due to the force. Fixing the bond 
length reduces the degrees of freedom to the orientations of the bonds, given 
by the polar angle 6'„ between the bond and the direction of the force and the 
azimuth angle (pn in a plane perpendicular to the force. So the hamiltonian 
reduces to 



6n,n+i is the angle between the bonds n and n + 1. We observe that, apart 
from a trivial additive term, the hamiltonian contains two dimensionless 
parameters 



m 




(86) 



y = biffu/T* 



and w = 6(/)//(2T*). 



(87) 



The partition function thus is obtained as 



- N-l N 

Z = Tr exp y ^ cos 9n,n+i + ^ cos 9. 



n 



(88) 



n=l n=l 
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where the trace is an integration over aU angles For the evaluation 

we use the expansion in terms of the Legendre polynomials 

eMyz) = lli'^i + ^)<ii{y)Pi{z), (89) 
I 

with the Pi{z) the Zth Legendre polynomial. The coefficients qi{y) are recov- 
ered as the integrals 

qiiy) = \ f_^dzeMy^) Pi{^)- (90) 

These elementary integrals, related to the modified half-integer Bessel-f unctions, 
can be calculated successively using the recurrence relations between the Leg- 
endre polynomials. A Legendre polynomial is subsequently expanded in the 
bond angles as 

{21 + 1) Pi{cOSen,n+l) = 47r ^y,,^(^„, 0„) Y,;„(^n+1, 4>n+l), (91) 

m 

with the Yi^ra as spherical harmonics. The advantage of this expression is 
that the relative angle between the bonds is expressed in terms of products 
of function of the the angles of the bonds in a fixed coordinate system. 

Next we formulate the partition function as a matrix product. The first 
step is the integration over the azimuth angles The angle 0i appears 
only in the first joint between 1 and 2. So its integration selects the term 
m = 0, as all the higher m integrate to zero. This implies that also the angle 
02 disappears from the product, since 

rz,o(^i,0i)>^*o(^2,02) = ^^PKcos^i)Pz(cos^2). (92) 
' 47r 

Subsequently one can integrate over 02 and so on, each integration over an 
azimuth angle contributing a factor 2tt. 

The remaining integration over the 9n is of the type with z = cos 6 

Tk,i{^w) = [(2fe + l)(2^ + l)]^/^ £ ^^p^^^^ exp{2wz) Pi{z). (93) 

The numerical factor in front of the integral is inserted to yield the property 

Tk,i {2w) = Y^ Tk,n {w) Tn,i {w) , (94) 

n 
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which follows from the completeness of the Legendre polynomials and which 
is convenient for symmetrizing the coming matrix product. 
We now define the transfer matrix 

Sk,i{w, y) = J2 Tk,n{w) Qniy) Tn,i{w), (95) 

n 

which allows us to write the partition function as 

Z=iAnf{x\S''-'\x), (96) 

with the definition of the initial state \x) 

\Xi) = Ti,o. (97) 

Note that we have used the property (194|) to make the definition of -B^,/ 
symmetric. 

By expression fl96p the problem of calculating the partition function is 
reduced to finding the eigenvalues of the matrix 5*^,;. For very long chains 
the largest eigenvalue will dominate. However, since the persistence length 
is long, this limit will only be reached when the chain has several thousands 
of bonds. Therefore we also investigate the finite chain for which one has to 
take into account all the eigenvalues. The eigenvalues and eigenfunctions are 
obtained from the equations 

E5mC = A«C. (98) 

n 

Then we use the representation 

5M = ECAaC. (99) 

a 

Note that for small w, i.e. for weak forces, the matrix Tkj reduces to 6k,i 
because the definition is just the orthogonality relation between Legendre 
polynomials. Then the eigenvalues of S^^i are given by the factor qi{y) and 
one of their properties is that they decrease with increasing /. Thus for w ~ 1, 
i.e. forces / ~ T*, it suffices to restrict the matrix 5* to a finite square in 
the upper-left corner. This simplifies the problem considerably since, both 
for the qi{y) and the Tkj{w), simple recursions relation exist, enabling to 
construct them from the first few values. Even for larger w one does not 
need to include an excessive number of states to get an accurate answer. 
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The partition function itself is not the main goal of the above exercise. 
We now apply this approach to the calculation of the force extension curve 
and the orientation correlation function. For the former we need the sum of 
the averages of cos6'j (see ( 120|) ). The appearance of cos6'j means the extra 
matrix tk,i at position j in the matrix product, with 

\(2k + l)(2l + , . X . ^ . X 

tk,i = '-^ L J_^dzPk{z)zPi{z). (100) 

This integral equals 

, _ + {k + l)4+i,f 

[(2A; + 1)(2/ + 1)]V2 • ^'^'^ 

As we plan to take the eigenfunctions S as basis, we transform the matrix 
tk^i to that basis: 

k,l 

In the same way we transform the x vector: 

k 

The partition function is now obtained as 

Z = (47r)^Ex"A^V. (104) 

a 

In the infinite long chain limit only the largest eigenvalues Aq contributes. 

The numerator of the end-to-end vector gives the sum over all positions 
of the cos6j, or matrices t in the product. We may write this as 



1 

a,l3 



(^cos%) = -5:x"f/.,/3t"'V. (105) 



The factor Ua,i3 accounts for all the ways the t can be put on the bonds. On 
the left of the t one will have the state a and on the right (3. So f/^ /3 is given 

by 

U.,, = X?-' + + ■■■ + A.A^^ + A^-\ (106) 

Dividing fllOSp by N will give the extension, i.e. the ratio between the end- 
to-end vector and the contour length. As mentioned in the limit of an infinite 
chain we may restrict the summations over a and /3 to a = /3 = 0. 
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For the orientation correlation we consider the infinitely long chain, as 
its definition should not be mudded by finite-size effects. We repeat the 
definition of the orientation correlation for the approximation considered in 
this appendix: 

Oij — (uj • Uj) — (cos 9i cos 9j + sin 9i sin 9j cos(0i — cos 0^)). (107) 

The first term gives the contribution of the longitudinal component, the sec- 
ond gives the transversal component (with respect to the force direction). 
We restrict ourselves to the first contribution, since we are mostly interested 
the persistence length of the force-free chain for which the eigenvalue prob- 
lem is trivial. Without an external force the longitudinal and transversal 
components are equal. The appearance of two terms cos 9i and cos 9j gives 
an insertion of the matrix tk,i at the bonds i and j. Outside the interval 
i — j the chain is in the state '0'^. Inside the interval it may be in any of the 
eigenstates of S. Thus we arrive at the expression for the average 

(cose, cos^,) = E/t°'"(Aa/Ao)l^-^lt'^'V- (108) 

a 

For the force-free chain this expression reduces considerably. As we men- 
tioned for w = the matrix T^^i reduces to 5k,i and each I is an eigenstate 
of Sk^i with eigenvalue qi{y). Therefore the t matrix with upper indices is 
the same as the one with lower indices. Moreover the t matrix couples only 
nearby states. So the summation over a is restricted to a = 1. The largest 
and the next eigenvalue read 

sinhi/ coshj/ sinhj/ 

Qoiy) = , qi{y) = —■ 109) 

y y y^ 

Therefore the expression for the orientation function reduces to 

(cose, cose,) = (gi(?/)Ao(z/))l'-^V3. (HO) 

One observes that the correlation decays exponentially with the persistence 
length Ip 



^P=- 1__/_ /,AN . /,A - (111) 



1 

log(?i(y)/?o(y)) log(cothy - 1/y) ' 

Large y gives the value Ip y. This has been used in the parameter match 
between our model and dsDNA. 
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— fit to experimental data, Wang et at. [2] 

— jointed chain model, N=°° 

— jointed chain model, N=64 
o simulation, our model, N=64 

— jointed chain model, N= 12S 
+ simulation, our model, = 128 

— jointed chain model, N=256 
X simulation, our model, N = 256 

5 10 15 20 

FU(kT) 

Figure 1: Force-extension curves for our model for = 64, 128 and 256, 
compared to the jointed chain model defined earlier in the main text for the 
same iV- values and to the experimental data by Wang et al. [7] . The agree- 
ment between our model (simulation) and the jointed chain (calculation) is 
excellent. Note that both models exhibit strong finite-size effects. 
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Figure 2: Eigenalues Aj, (dotted lines) and A* (solid lines and points) as 
a function of mode number p. Here, the model parameters u = 0.35 and 
T* = 0.034 are chosen to represent dsDNA. The polymer lengths are, from 
top to bottom, = 64, 128 and 256. The black lines correspond to Eq. (1811) 
for the respective N values. 
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